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Abstract. We introduce a simple model which shows non-trivial self organized critical properties. The 
model describes a system of interacting units, modelled by Polya urns, subject to perturbations and which 
occasionally break down. Three equivalent formulations - stochastic, quenched and deterministic - are 
shown to reproduce the same dynamics. Among the novel features of the model are a non-homogeneous 
stationary state, the presence of a non-stationary critical phase and non-trivial exponents even in mean 
field. We discuss simple interpretations in term of biological evolution and earthquake dynamics and we 
report on extensive numerical simulations in dimensions d — 1, 2 as well as in the random neighbors limit. 

PACS. 64.60.Ht Dynamic critical phenomena - 64.60.Lx Self-organized criticality; avalanche effect 



Our understanding of Self Organized Criticality (SOC) 
|0| , as a general framework for the emergence of scale- free 
behavior in Nature, has greatly benefitted from the in- 
troduction of simple models. Even though models such 
as the sandpile J5J and the Bak-Sneppen [B are too sim- 
ple to capture the complexity of natural phenomena such 
as earthquakes || and biological evolutional^], they have, 
nonetheless, identified some basic mechanisms leading to 
SOC. These systems have been a starting point both for 
the development of more complex and realistic models of 
natural phenomena and for analytical approaches |^|)]. 
which have led us to a much deeper understanding of SOC. 
Indeed, we can now identify some basic "routes to Self Or 
■ ganized Criticality" such as those based on sandpile [p] 
\ extremal dynamics |l0| , |r^| , memory ||l3f and network 
models. 

In this Rapid Communication we propose a qualita- 
tively different "route to SOC" based on a very simple 
model of interacting Polya urns. Its qualitative differences 
with respect to other SOC models are that it is character- 
ized by a non-homogeneous stationary state and by non 
trivial exponents even in the mean field case. Furthermore, 
we shall show numerical evidence for the occurrence of a 
non- stationary self organized critical state. Moreover, the 
model can be formulated in three different but equivalent 
ways. This fact, on one hand allows us to use a wide vari- 
ety of tools to investigate its critical properties, and on the 
other it bridges different descriptions of the same process. 
All these features can well be relevant in the description of 
natural phenomena. The model indeed provides a general 
framework for the emergence of SOC which, as we shall 
discuss, can be applied both to coevolution and to large 
scale earthquakes dynamics. Note indeed that the pat- 
terns of earthquakes activity are highly non-homogeneous 



and that such a system is, in principle, non-stationary. 
The same applies to our ecosystem, which is in a non- 
stationary state where ever fitter species replace less fit 
ones. 

We consider a system of interacting Polya urns ar- 
ranged on a d-dimensional lattice. A Polya urn is a simple 
model to study e.g. the occurrence of accidents |Q. Each 
urn contains initially b black balls and 1 white one. As 
_p_sandpile models, at each time step we randomly select 
ilOlite and attempt to add a "grain of sand" , i.e. a white 
ball, to the corresponding urn. A ball is drawn from the 
selected urn: If the ball is white the attempt is successful 
and a new white ball is added to the urn. If it is black 
a "fatal accident" occurs: The urn becomes unstable and 
it "topples". The toppling mechanism is as follows: 1) the 
urn is reset to 1 white ball and b black ones and 2) for each 
white ball of the original urn a similar attempt to add a 
white ball is made on a randomly chosen nearest neighbor 
urn. In this way, white balls released by an unstable urn 
can provoke some "fatal accident" in nearby urns (addi- 
tion of white ball to already unstable urns leaves them 
unstable but it increases the number of white balls in it). 
The process stops when all balls are redistributed provok- 
ing no further toppling. A new attempt to add a ball to a 
randomly chosen urn is made, at the next time-step, and 
the process goes on. Dissipation of balls at the boundary, 
as in the sandpile Q, can also be considered to allow the 
system to relax to a stationary state. Actually, in order 
to keep the same definition of the model both in finite 
dimensions and in the random neighbor version, we con- 
sider here "bulk dissipation" modifying step 2) into: 2') 
with probability A all white balls disappear, otherwise 2) 
applies. 
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Fig. 1. Snapshot of a section of the d 
L = 128 for b = 2.50 and A = 10" 2 . 



2 system of size 



Let e, be the number of subsequent additions of white 
balls in urn i, since the last draw of a black one. Then urn 
i contains b black balls and ej + 1 white ones. The toppling 
probability, i.e. the probability to draw a black ball from 
urn i, is then 

e % + b + 1 

In general, the probability that urn i topples, if it receives 
q white balls from neighbors, is 



Wi(q) 



1 



(ej+qy.rje. + b+l) 
e^.r^ + q + b + l) 



The -T-function is introduced here in order to generalize 
our discussion to any real positive b. 

The system spontaneously evolves to a critical state 
which is generally characterized by a non uniform distri- 
bution of the variables ej. A snapshot of the system for 
b = 2.5 is shown in Fig. [I]. It is clearly seen that very 
stable urns 3> 1 coexist with less stable ones. 

We say that the initial perturbation causes an avalanche 
of size s and volume v, where s is the number of toppling 
events occurring before the system returns stable and v is 
the number of distinct sites involved in the avalanche. Op- 
erationally all urns which become unstable after the first 
event topple simultaneously. The urns, which as a result 
of this first wave of topplings become unstable, topple si- 
multaneously in a second wave, and so on. An avalanche is 
then also characterized by the number w of such waves oc- 
curring before the system returns stable. Finally, one can 
also measure the total number E of white balls involved 
in the avalanche event. Of course, if the initial urn resists 
absorbing an extra white ball, s = v = w = E = 0. The 
distribution of avalanche sizes in the critical state 



has a power law behavior. In this state the volume, waves 
and the number E of the avalanche are related to its size 
by power laws 



E 



(4) 



P(s) 



(3) 



The four exponents r, z, u> and 7 define the critical state 
of the model. They are given in Table |l| for several values 
of b > 1 in d = 1, d = 2 and in the random neighbors 
version (K "neighbors" are randomly chosen each time 
among all sites) respectively^} For b < 1 the system be- 
comes non-stationary: The number of white balls in the 
system increases linearly in time. In spite of this, we have 
found that avalanches have still a well defined distribu- 
tion. We have analyzed in particular the border-line case 
6=1: For the random neighbor model, we have found that 
in a range of times (from 10 5 to 2 • 10 8 ), for various system 
sizes (up to n — 2 14 urns) and different dissipation rates 
(A = 1CP 2 and A = 10 -3 ), the distribution of avalanche 
sizes follows Eq. (^) on more than three decades, with an 
exponent which is, for all these cases, always in the range 
r G [1.95,2.05]. In d = 2, for sizes up to n = 2 14 and in 
a range of times from 5 ■ 10 7 to 3 ■ 10 s , we have similarly 
found r £ [2.03, 2.11]. In d = 1, extensive simulations over 
lattice lengths of 256 and 512 sites and on a range of times 
from 10 9 to 6 • 10 9 , we have found t £ [2.10, 2.21]. 

In an ecosystem species are probed by changes in the 
environment in a fashion which goes under the name of co- 
evolution. Namely, the environment, which is constituted 
by the interaction among all living organisms, is implicitly 
modified by each single being and determines its fate as 
well^j. From the point of view of any single species, such 
perturbations impose a random selective pressure which 
eventually leads to a change in the constitution of the 
species or to its extinction. A particular character, devel- 
oped by a random mutation in a subspecies, can be "se- 
lected" by evolution if it becomes essential for the survival 
of the whole species in a changed environment j|] . In this 
process, which is driven by chance, species become more 
and more complex. A more complex species has also a 
wider variability in subspecies which, in this process, will 
make it more resistant. In much the same way, Polya urns 
in the model "evolve" increasing their "complexity" and 
their robustness to external perturbations. The model also 
assumes that the extinction of well adapted "species" (i.e. 
urns with > 1) produces a larger perturbation than 
that caused by poorly fit species. 

In this toy ecosystem, each species has exactly the 
same chances of any other species to survive, when it ap- 
pears (e^ = 0). There is no a priori genetic character- 
istic, such as fitness, which guarantees the survival of a 

1 The parameter A, which introduces bulk dissipation, sets 
the upper cutoff of the critical states. We found that, in d = 
1,2, dissipation at the boundaries, as in the sandpile model, 
leads to similar results. In the random neighbors version, we 
have found that the exponents do not depend on the number 
K of neighbors chosen. 

2 Just as an example, one could think at the selection pres- 
sure produced by the precence of oxygen in the air, which has 
been initially produced by some zoo-phyte. 
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Table 1. The exponents for different values of b for d — 1 and sizes up to L = 512, d = 2 and L = 128, and for the random 
neighbor model with n — 16384 sites and K = 2. 



species. Its survival will rather depend on its "ability" to 
adapt constantly, via random mutations, to the changing 
environment. This perspective, also suggests that an op- 
erational measure of fitness (or resistance) of a species, is 
possible using Eq. (|l|): High ej means high fitness. Note 
that this differs from the concept of fitness introduced in 
most SOC models of coevolution[p[|^,|l4| . In these, fitness 
is related to reproduction rates, whereas in our simplified 
picture of coevolution, fitness emerges as a measure of the 
resistance against extinction of a species. 

We show now that this different notion of fitness, when 
introduced as an intrinsic property of each species, leads 
to exaclty the same coevolutionary process. To be more 
precise, let us assume that the probability /$ of extinction 
of species i under a perturbation, is no more given by Eq. 
(|l|), but it is rather fixed for each species. In particular, 
this intrinsic property is drawn randomly for each species 
from a distribution p(f). Accordingly we also replace Eq. 
(||) by Wi(q) = 1 — (1 — fi) q . Species are probed by ran- 
dom perturbations (addition of white balls) and, as before, 
the extinction of species i perturbs "neighboring" species 
in the interaction web via the same toppling mechanism. 
When a species disappears, its "niche" (site) is immedi- 
ately occupied by a new species, with a new randomly 
drawn fitness value f[. Thus, for 



b-l 



(5) 



we obtain exactly the same stochastic dynamics given by 
Eqs. (0,0). In order to show the equivalence, it is enough 
to show that Eq. (p|) leads to the same rates Wi (q) of Eq. 
(|J). Consider one particular urn and let e^t be the value 
of the corresponding variable after t drawings. The event 
£i,t — e occurs with a probability ^ P(ei.t — e\ei,t-e = 
0, fi) = (1 — fi) e . Taking the average over p(f), we find 



ee«,t- 



0) 



r(6 + l)e! 
r(e + b+l) 



(6) 



Using detailed balance, P(ei t t+ q = e + q\ei tt - e = 0) = 
Wi(q)P(ei lt = e\e iit -e = 0), we easily recover Eqs. (Q, §). 

3 Here t is a "local" time on site i, which measures the num- 
ber of perturbations on that site. Therefore the event a t t = e 
implies ei t t- e = 0. Accordingly we used the notation P(ei t t = 
e\e it t- e = 0,/i) (note indeed that P(e itt = 0|e <jt = 0, fi) = 1 
for e = 0). 



The equivalence of the two models has been also tested in 
numerical simulations. 

The initial definition of the model is completely stochas- 
tic, whereas in the formulation based on Eq. (g|), fi are 
fixed, quenched variables, which are renewed stochasti- 
cally at each extinction event. The equivalence of the two 
descriptions is an example of a general mapping p"2|, re- 
cently developed to deal with extremal dynamics. Its ap- 
plication in the biological context are also discussed in 
Ref. @. 

There is a further interesting mapping, originally de- 
veloped in the context of interface growth[[l^, which can 
be applied to the present model. As in the sandpile modcljlj 
we define the toppling probability as 



fi = if ej < hi and fi — I if e, > hi 



(7) 



While in the sandpile the thresholds are fixed hi = 2d, 
we introduce here a model where hi are randomly drawn 
from a given distribution ip(h) after each toppling event 
on site i. The choice 



m = 



br(b + i)r(h) 
r(h + b + i) 



(8) 



reproduces a dynamics which is equivalent to the previous 
two formulations of the model. To prove this, it is enough 
to derive the statistics of the number ej = hi — 1 of per- 
turbations that an urn with fi = f will overcome before 
toppling. Clearly P{h % = h\fi= f) = f(l - f) h ~\ Taking 
the average over the distribution (||) of / leads indeed to 
Eq. ®. 

This formulation is completely deterministic: It as- 
sumes that each urn appears with a prescribed "lifetime" 
measured in terms of perturbations. As soon as this life- 
time is reached, the urn topples. This is the same thresh- 
old dynamics used in the sandpile model. Here however 
thresholds hi are very broadly distributed (note that ^j(h) ^ 
/i _b_1 ), whereas in the sandpile ip(h) = 5(h — 2d). The 
sandpile is a paradigm for seismic phenomena: Each site 
represent a fault which is perturbed by the slow addition 
of stress energy. When the energy load of a site exceeds 
the threshold hi, the fault breaks down and all the energy 
is released to neighbor sites. Our model also proposes a 
different description of the same phenomenon: each addi- 
tion of stress energy has the same probability /j to provoke 
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an earthquake. When this occurs, it provokes the energy 
release and a local seismic rearrangement, i.e. fi — > f[. 

Using the quenched version of the model given by Eq. 
(||), it is easy to derive the effective distribution p(f) of 
the fi's in the system at the stationary state. It is enough 
to consider detailed balance in an interval fi G [/, / + 
df) under a single perturbation. The probability that, in 
one time step, one fi leaves this interval is fp(f)df. This 
has to balance the number p(f)df f Q df'f'p(f') of sites 
which enter this interval. This gives p(f) = (b — l)f b ~ 2 . 
With Eq. (H) one can also derive the distribution of 6j 
in the system. Indeed, P(ei. t = e) = P(e-i,t — e\ei,t-e = 
0)P(e M _ e = 0) where P(e t = 0) = 1 - 'l/b is derived 
imposing normalization. This leads to -P(ej) ~ e~ b . 

As b — > 1 + , both the distributions p(f) and P{ei) be- 
come unnormalizable and the probability of finding sites 
with a = vanishes. Accordingly, numerical simulations 
show that, for b < 1, the system average of &i increases 
linearly with time and the system never reaches a sta- 
tionary state. The divergence of normalization of p(f) at 
/ = occurs because less "fit" species are more rapidly 
replaced than more "fit" ones. For b < 1 the search for 
the "perfect" species fi = never stops. For b > 1, 
the probability that a perturbation causes a toppling is 
P top = J dffp(f) = which also vanishes as b — ► 1. 
This means that, for b < 1, avalanches occur more and 
more rarely as time goes on. 

Let us discuss in more detail the random neighbor 
model. Numerical results are consistent with z/D = 1 and 
7 = max[l, l/(b— 1)], which is what one expects from the 
observation that each site is involved at most once in the 
same avalanche (the exponent 7 = 1/(6 — 1) for 1 < b < 2 
comes from the limit laws of Levy variables). On the other 
hand we see that the exponents r and to differ from their 
usual mean-field values r = 3/2 and u> — 1/2, and that 
such values are eventually reached for b — -> 00. This devia- 
tion r ^ 3/2 can be understood as an effect of correlation. 
In order to show this, let us review the argument leading 
to r — 3/2. Consider an avalanche and let M t be the num- 
ber of unstable sites after t toppling events. M t performs 
a random walk and the size s = min{£ : M t = 0, t > 0} 
of the avalanche is the first return time of M t to 0. There- 
fore s has the same distribution of the first return times 
to the origin of a random walk P(s) ~ s -3 / 2 . Since the 
steps \M t — Mt-x\ of the random walk are bounded by 
the coordination number K, the only possibility for a de- 
viation from r = 3/2 is to have correlations. Correlations 
indeed arise because a toppling event may release many 
white balls and these are transported along the avalanche 
thus modifying the probability of toppling of successive 
sites. The theoretical calculation of the exponents for the 
random neighbor version is a challenging problem under 
current investigation. 

In conclusion, we have introduced a very simple model 
which displays non-trivial self organized critical features. 
The model was analyzed numerically and we also derived 
some analytic result. The main features of the model are 
non-homogeneous critical states, a critical non-stationary 
state in a region of the control parameter (b < 1) and 



non-trivial exponents even in the mean field limit. Further- 
more, it allows for three different equivalent formulations, 
which allow one to better investigate and understand the 
nature of the critical state. As a closing remark, we notice 
that the non-stationary regime can be relevant both for 
the description of ecological systems and for earthquake 
dynamics. Both these systems are indeed not in a station- 
ary state. Interestingly enough, the exponent r is larger 
than 2 for b < 1. Both the Gutenberg- Richter law for 
earthquakes H and the extinction size distribution)^] also 
display avalanche exponents close to 2. 
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